#Replication File for
#Wang, Austin Horng-En. 2025. 
#"Political Opportunity Structure Conditions the Legacy of Political Violence",
#Journal of Conflict Resolution, Accepted
#Version: March 30, 2025

#Study 2: Individual-level Evidence

rm(list=ls(all=TRUE))

#Create summarySE function for summary
summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
                      conf.interval=.95, .drop=TRUE) {
  library(plyr)
  length2 <- function (x, na.rm=FALSE) {
    if (na.rm) sum(!is.na(x))
    else       length(x)
  }
  datac <- ddply(data, groupvars, .drop=.drop,
                 .fun = function(xx, col) {
                   c(N    = length2(xx[[col]], na.rm=na.rm),
                     mean = mean   (xx[[col]], na.rm=na.rm),
                     sd   = sd     (xx[[col]], na.rm=na.rm)
                   )
                 },
                 measurevar
  )
  
  datac <- rename(datac, c("mean" = measurevar))
  
  datac$se <- datac$sd / sqrt(datac$N)  # Calculate standard error of the mean
  
  ciMult <- qt(conf.interval/2 + .5, datac$N-1)
  datac$ci <- datac$se * ciMult
  
  return(datac)
}

#Call data in. Please replace the link to the your local folder
library(foreign)
library(car)
project1 = read.spss("C:/Users/weary/Dropbox/Chen_Wang_Wu_Yeh_Projects/referendum proposal/PP2197EC2/PP2197EC2/PP2197EC2.sav",
                     to.data.frame=TRUE,use.value.labels=FALSE)

#Votechoie for all respondents

project1$votechoice=NA
project1$votechoice[project1$Q24a==1]="KMT"
project1$votechoice[project1$Q24b==1]="KMT"
project1$votechoice[project1$Q24c==1]="KMT"
project1$votechoice[project1$Q24a==2]="DPP"
project1$votechoice[project1$Q24b==2]="DPP"
project1$votechoice[project1$Q24c==2]="DPP"


table(project1$votechoice)
#445DPP 447KMT

#Three Treatment (Tie, KMT55, DPP55)

project1$group3=NA

project1$group3[project1$Q24a<=2]="Tie"
project1$group3[project1$Q24b<=2]="KMT55"
project1$group3[project1$Q24c<=2]="DPP55"

table(project1$group3)  
#DPP55 KMT55   Tie    #The number of assignments in p.25
#305   299   298 


#Anova test for votechoice and the treatment (page 26)
table(project1$group3,project1$votechoice)
prop.table(table(project1$group3,project1$votechoice),1)
chisq.test(table(project1$group3,project1$votechoice))
#p=0.044 (p.26)

#The White Terror Priming  (p.22)
table(project1$Q22)
project1$treatment2=NA

project1$treatment2[project1$Q22<91]="Treatment"
project1$treatment2[project1$Q22==99]="No T"
table(project1$treatment2)

#No T Treatment    (p. 23, half of respondents received the treatment)
#456       454 

table(project1$treatment2,project1$votechoice)

#Examine the interaction effect, p.28
"the KMT-dominance in the pre-election poll treatment significantly increased 
#the likelihood for Taiwanese respondents to support KMT (p = 0.036. 
#KMT received 48.6%, 58.7%, and 44.1%, respectively), "

table(project1$group3[project1$treatment2=="Treatment"],
      project1$votechoice[project1$treatment2=="Treatment"])

prop.table(table(project1$group3[project1$treatment2=="Treatment"],
                 project1$votechoice[project1$treatment2=="Treatment"]),1)

chisq.test(table(project1$group3[project1$treatment2=="Treatment"],
                 project1$votechoice[project1$treatment2=="Treatment"]),1)
#p-value = 0.03586 (p.28)

#No difference for no white terror treatment group (p.27)
#One-way ANOVA shows no significant difference across the three groups 
#on the left (p = 0.57. KMT received 49.7%, 51.3%, and 45.3%, respectively). 
table(project1$group3[project1$treatment2=="No T"],
      project1$votechoice[project1$treatment2=="No T"])

prop.table(table(project1$group3[project1$treatment2=="No T"],
                 project1$votechoice[project1$treatment2=="No T"]),1)

chisq.test(table(project1$group3[project1$treatment2=="No T"],
                 project1$votechoice[project1$treatment2=="No T"]))
#p = 0.5701 (p.27)


#=================================================#
#Figure 4  (p.27)

project1$votechoice2=NA
project1$votechoice2[project1$votechoice=="KMT"]=1
project1$votechoice2[project1$votechoice=="DPP"]=0
table(project1$votechoice2)

result1=summarySE(project1,"votechoice2",groupvars=c("treatment2","group3"))
result1=result1[complete.cases(result1),]
result1

library(ggplot2)

result1$treatment2[result1$treatment2=="Treatment"]="White Terror\nTreatment"
result1$treatment2[result1$treatment2=="No T"]="No Treatment"


ggplot(data=result1,aes(x=treatment2,y=votechoice2,fill=group3))+
  geom_bar(position=position_dodge(.9), stat="identity", colour='black') +
  geom_errorbar(aes(ymin=votechoice2-ci, ymax=votechoice2+ci), width=.2,position=position_dodge(.9))+
  scale_fill_manual(labels=c("DPP55","KMT55","No Info"),values=c("green2","blue","gray"))+theme_bw()+
  guides(fill=guide_legend(title="Exp Group"))+ylim(0,1)+ylab("% support KMT")+xlab("")


ggplot(data=result1,aes(x=treatment2,y=votechoice2,fill=group3))+
  geom_bar(position=position_dodge(.9), stat="identity", colour='black') +
  geom_errorbar(aes(ymin=votechoice2-ci, ymax=votechoice2+ci), width=.2,position=position_dodge(.9))+
  scale_fill_manual(labels=c("DPP55","KMT55","No info"),values=c("green2","blue","gray"))+theme_bw()+
  guides(fill=guide_legend(title="Exp Group"))+ylim(0,1)+ylab("% support KMT")+xlab("")+
  geom_segment(aes(x=0.7, xend=1.3, y=0.8, yend=0.8),colour="black") + 
  geom_segment(aes(x=1.7, xend=2.3, y=0.8, yend=0.8),colour="black") +
  geom_text(x=1,y=0.83,label="ANOVA, p = 0.57",colour="black",size=3.5)+
  geom_text(x=2,y=0.83,label="ANOVA, p = 0.036*",colour="black",size=3.5)

ggsave("Figure 4.png",width= 6, height = 4, units="in",dpi=300)

#=============================================#
#Randomization check and Control Variales (p.26)
#And Table 2 (p.23)
project1$group6=NA
project1$group6=paste(project1$treatment2,project1$group3)
table(project1$group6)
project1$group6[project1$group6=="No T NA"]=NA
project1$group6[project1$group6=="Treatment NA"]=NA

table(project1$group6)

#Female
table(project1$sex)
project1$female=NA
project1$female=recode(project1$sex,"1=0;2=1")
#age
table(project1$age)
#edu
table(project1$edu)
#Taiwan Identity
table(project1$Q26)
project1$TWID=recode(project1$Q26,"1=1;else=0")
#Partisanship
table(project1$Q27)
project1$BLUE=recode(project1$Q27,"1=1;3=1;5=1;6=1;else=0")
project1$GREEN=recode(project1$Q27,"2=1;4=1;7=1;8=1;9=1;else=0")
#nonepartisans
project1$NONE=1
project1$NONE[project1$BLUE==1]=0
project1$NONE[project1$GREEN==1]=0
table(project1$NONE)
summary(project1)

#Randomization check, by treatments
#One-way ANOVA shows no significant difference for the two treatments 
#in the respondents�� gender  (p.26)
summary(aov(sex~treatment2,data=project1)) #p=0.881
summary(aov(sex~group3,data=project1)) #p=0.297

#age
summary(aov(age~treatment2,data=project1))  #p=0.329
summary(aov(age~group3,data=project1)) #p=0.928

#education
summary(aov(edu~treatment2,data=project1))  #p=0.196
summary(aov(edu~group3,data=project1))  #p=0.687

#TWID
summary(aov(TWID~treatment2,data=project1))  #p=0.351
summary(aov(TWID~group3,data=project1))  #p=0.259

#Two way anova (p.26)
summary(aov(age~group3*treatment2,data=project1)) #p=0.414 for age
summary(aov(edu~group3*treatment2,data=project1)) #p=0.254 for edu
summary(aov(TWID~group3*treatment2,data=project1))#p=0.216 for TWID 
summary(aov(female~group3*treatment2,data=project1)) #p=0.442 for female


#=============================================#
#Regression Analysis, Table 3, p. 28
project1$group3=relevel(as.factor(project1$group3),ref="Tie")


model1=glm(votechoice2~group3,data=project1[project1$treatment2=="Treatment",],family='binomial')
summary(model1)

model2=glm(votechoice2~group3+female+TWID+age+edu,data=project1[project1$treatment2=="Treatment",],family='binomial')
summary(model2)

model3=glm(votechoice2~group3,data=project1[project1$treatment2=="No T",],family='binomial')
summary(model3)

model4=glm(votechoice2~group3+female+TWID+age+edu,data=project1[project1$treatment2=="No T",],family='binomial')
summary(model4)

library(stargazer)  #Table 3
stargazer(model1,model2,model3,model4,type="text")


#============================================#
#Final robustness check (p.29)

chisq.test(table(project1$treatment2[project1$group3=="KMT55"],
                 project1$votechoice[project1$group3=="KMT55"]))
#p=0.25, chi square test, for two KMT55 groups

chisq.test(table(project1$treatment2[project1$group3=="KMT55"],
                 project1$votechoice[project1$group3=="KMT55"]))

table(project1$votechoice2,project1$group6)
#KMT55+white terror treatment group is statistically different from the two control groups
chisq.test(matrix(c(60, 85, 75, 62), nrow = 2)) #p=0.033 for control + no treatment
chisq.test(matrix(c(60, 85, 90, 71), nrow = 2)) #p=0.11 for control +treatment
chisq.test(matrix(c(60, 85, 73, 69), nrow = 2)) #p=0.02 for control +treatment

           